#!/usr/bin
# -*- coding: utf-8 -*-

from math import sqrt as msqrt

EPSILON = 1.0e-6
ABAQUS_MAX_ARGS_PER_LINE = 16

fileName = "maillage2.dat"

nodes = {}
elements = {}
domains = {}
topologicDomains = {}

# Lecture du fichier d'entrée ==================================================
try:
    myFile = open(fileName, "r")
except:
    print "Can not open " + fileName

iLine = 0
while True:
    line = myFile.readline()
    if not line:
        break

    iLine += 1
    line = line.replace("\n","")
    line = line.replace("\r","")
    line = line.split(",")

    if line[0].lower() == "n":
        try:
            nId, nx, ny, nz = line[1:]
        except:
            raise ValueError, "Line " + str(iLine) + ": reading node failure"
        try:
            a = nodes[nId]
        except:
            nodes[nId] = (float(nx), float(ny), float(nz))

    elif line[0].lower() == "e":
        try:
            eSection, eId, eNds = line[1], line[2], line[3:]
        except:
            raise ValueError, "Line " + str(iLine) + ": reading element failure"
        try:
            a = elements[eId]
        except:
            elements[eId] = [eSection]
            for n in eNds:
                elements[eId].append(n)

    elif line[0].lower() == "d":
        try:
            dId, de = line[1], line[2:]
        except:
            raise ValueError, "Line " + str(iLine) + ": reading domain failure"
        try:
            a = domains[dId]
            for e in de:
                domains[dId].append(e)
        except:
            domains[dId] = de

    elif line[0].lower() == "td":
        try:
            tdId, tdn = line[1], line[2:]
        except:
            raise ValueError, "Line " + str(iLine) + ": reading domain failure"
        try:
            a = topologicDomains[tdId]
            for n in tdn:
                topologicDomains[tdId].append(n)
        except:
            topologicDomains[tdId] = tdn

    elif line[0] == "":
        pass

    else:
        raise ValueError, "Line " + str(iLine) + " must start by: n, e or d"

myFile.close()

print "Number of nodes............: " + str(len(nodes))
print "Number of elements.........: " + str(len(elements))
print "Number of domains..........: " + str(len(domains))
print "Number of topologic domains: " + str(len(topologicDomains))

if False:
    for key, val in nodes.items():
        print "node", key, val

    for key, val in elements.items():
        print "element", key, val


# Connectivité complémentaire ==================================================
nodesFromDomain = {}
for dId, dEls in domains.items():
    nodesFromDomain[dId] = {}
    for eId in dEls:
        nds = elements[eId]
        for nId in nds:
            try:
                a = nodesFromDomain[dId][nId]
            except:
                nodesFromDomain[dId][nId] = None
    nodesFromDomain[dId] = nodesFromDomain[dId].keys()



# Connectivité inverse =========================================================
elementsFromNode = {}
domainsFromElement = {}
domainsFromNode = {}

for elId, elVal in elements.items():
    elSection = elVal[0]
    elNds = elVal[1:]

    for nd in elNds:
        try:
            a = elementsFromNode[nd]
            try:
                b = a[elId]
            except:
                elementsFromNode[nd][elId] = None
        except:
            elementsFromNode[nd] = {elId:None}


for dId, dEls in domains.items():
    for elId in dEls:
        try:
            a = domainsFromElement[elId]
            try:
                b = a[dId]
            except:
                domainsFromElement[elId][dId] = None
        except:
            domainsFromElement[elId] = {dId:None} # pour l'avenir avec des éléments pouvant appartenir à plusieurs domaines

        for nId in elements[elId][1:]:
            try:
                a = domainsFromNode[nId]
                try:
                    b = a[dId]
                except:
                    domainsFromNode[nId][dId] = None
            except:
                domainsFromNode[nId] = {dId:None}

for key, val in topologicDomains.items():
    print "Topologic domain:", key, "; corner nodes:", val


# Frontières des domaines topologiques =========================================
topologicEdges = {}
for tdId, tdCornerNodes in topologicDomains.items():
    topologicEdges[tdId] = []
    for i, firstEdgeNode in enumerate(tdCornerNodes):
        try:
            a = [int(firstEdgeNode),int(tdCornerNodes[i+1])]
            a.sort()
            a = [str(a[0]), str(a[1])]
            topologicEdges[tdId].append(tuple(a))
        except:
            a = [int(firstEdgeNode),int(tdCornerNodes[0])]
            a.sort()
            a = [str(a[0]), str(a[1])]
            topologicEdges[tdId].append(tuple(a))
            topologicEdges[tdId].append(tuple(a))

for key, val in topologicEdges.items():
    print "Topologic domain:", key, "; edges:", val


topologicEdgeNodes = {}
for tdId, tdEdgeCorners in topologicEdges.items():
    topologicEdgeNodes[tdId] = {}
    for ecId in tdEdgeCorners:
        n1, n2 = ecId
        # Vecteur de l'arête
        x1, x2 = nodes[n1], nodes[n2]
        u12 = (x2[0]-x1[0], x2[1]-x1[1], x2[2]-x1[2])
        nu12 = msqrt(u12[0]**2 + u12[1]**2 + u12[2]**2)

        topologicEdgeNodes[tdId][ecId] = {n1:0.0, n2:nu12}

        treatedNodes = {n1:None, n2:None}
        for nId in nodesFromDomain[tdId]: # attention, il faut que les identifiants coïncident entre domain et topologicDomain
            try:
                a = treatedNodes[nId]
            except:
                xn = nodes[nId]
                u1n = (xn[0]-x1[0], xn[1]-x1[1], xn[2]-x1[2])

                # Produit vectoriel pour savoir si les vecteurs sont colinéaires
                u12vecu1n = (u12[1]*u1n[2] - u12[2]*u1n[1],
                             u1n[0]*u12[2] - u1n[2]*u12[0],
                             u12[0]*u1n[1] - u12[1]*u1n[0])
                nu12vecu1n = msqrt(u12vecu1n[0]**2 + u12vecu1n[1]**2 + u12vecu1n[2]**2)

                if nu12vecu1n > EPSILON:
                    # Vecteur non nul ==> il n'appartient à la droite
                    pass
                else:
                    # Vecteur nul ==> il appartient à la droite il faut vérifier qu'il appartient au segment
                    u12scau1n = u12[0]*u1n[0] + u12[1]*u1n[1] + u12[2]*u1n[2]
                    if u12scau1n > 0:
                        # il est du bon côté du premier noeud
                        u2n = (xn[0]-x2[0], xn[1]-x2[1], xn[2]-x2[2])
                        u12scau2n = u12[0]*u2n[0] + u12[1]*u2n[1] + u12[2]*u2n[2]
                        if u12scau2n < 0:
                            # le noeud est bon
                            nu1n = msqrt(u1n[0]**2 + u1n[1]**2 + u1n[2]**2)
                            topologicEdgeNodes[tdId][ecId][nId] = nu1n
                treatedNodes[nId] = None

for key1, val1 in topologicEdgeNodes.items():
    print "topologicDomain:", key1
    for key2, val2 in val1.items():
        print "- arête:", key2
        for key3, val3 in val2.items():
            print "  - noeud:", key3, " ; abscisse curviligne:", val3


boundaries = {}
for td1, edges1 in topologicEdgeNodes.items():
    for td2, edges2 in topologicEdgeNodes.items():
        if td1 != td2:
            for edge1, nodes1 in edges1.items():
                for edge2, nodes2 in edges2.items():
                    if edge1 == edge2:
                        try:
                            a = boundaries[(td1,td2)]
                        except:
                            boundaries[(td1,td2)] = {edge1:(nodes1,nodes2)}

for key, val in boundaries.items():
    print "boundary:", key,
    print "- edge:", val.keys()
    print "- nodes:", val.values()




# Relations interdomaines ======================================================
domainsConnections = {}
for d in domains.keys():
    domainsConnections[d] = {}


# Boucle sur les noeuds
for nId, dIds in domainsFromNode.items():
    dIds = dIds.keys()

    # Si le noeud est lié à plusieurs domaines, alors il est un noeud frontière
    if len(dIds) > 1:

        # Boucle sur les domaines partageant le noeud
        for i, dId in enumerate(dIds):

            # Boucle sur les domaines liés au domaine d'identifiant dId
            otherdIds = dIds[:i] + dIds[i+1:]
            for otherdId in otherdIds:
                try:
                    c = domainsConnections[dId][otherdId]
                    try:
                        d = c[nId]
                    except:
                        domainsConnections[dId][otherdId][nId] = None

                except:
                    domainsConnections[dId][otherdId] = {nId:None}

if False:
    for key1, val1 in domainsConnections.items():
        print "Domain " + key1 + " is connected to domains:"
        for key2, val2 in val1.items():
            print "- " + key2 + " via nodes:"
            print "  - " + str(val2.keys())
        print "\n" * 3



# Génération d'un fichier inp ==================================================
myFile = open("maillage2.inp", "w")

# Header
myFile.write("*Heading\n")
myFile.write("*Part, name=plaque\n")
myFile.write("\n" * 2)

# Noeuds
myFile.write("*node, nset=allNodes\n")
for nId, nX in nodes.items():
    myFile.write("  " + nId + "," + str(nX[0]) + "," + str(nX[1]) + "," + str(nX[2]) + "\n")
myFile.write("\n" * 2)

# Éléments
myFile.write("*element, type=S4, elset=allElements\n")
for eId, eArgs in elements.items():
#    myFile.write("  " + eId + "," + eArgs[1] + "," + eArgs[2] + "," + eArgs[3] + "," + eArgs[4] \
#                            + "," + eArgs[5] + "," + eArgs[6] + "," + eArgs[7] + "," + eArgs[8] \
#                            + "\n")
    myFile.write("  " + eId + "," + eArgs[1] + "," + eArgs[2] + "," + eArgs[3] + "," + eArgs[4] \
                            + "\n")
myFile.write("*End Part\n")
myFile.write("\n" * 2)

# Assemblage
myFile.write("*Assembly, name=Assembly\n")
myFile.write("*Instance, name=plaque-1, part=plaque\n")
myFile.write("*End Instance\n")
myFile.write("\n" * 2)

# Domaines
# - éléments
for dId, dArgs in domains.items():
    myFile.write("*Elset, elset=domain_" + dId + ", instance=plaque-1\n")

    n = len(dArgs) // ABAQUS_MAX_ARGS_PER_LINE
    r = len(dArgs) % ABAQUS_MAX_ARGS_PER_LINE

    for i in xrange(n):
        line = "  "
        for j in xrange(ABAQUS_MAX_ARGS_PER_LINE):
            line += str(dArgs.pop()) + ","
        line = line[:-1] + "\n"
        myFile.write(line)
    line = " "
    for i in xrange(r):
        line += str(dArgs.pop()) + ","
    line = line[:-1] + "\n\n"
    myFile.write(line)

# - noeuds
for dId, dNds in nodesFromDomain.items():
    myFile.write("*NSet, nset=domain_" + dId + ", instance=plaque-1\n")

    n = len(dNds) // ABAQUS_MAX_ARGS_PER_LINE
    r = len(dNds) % ABAQUS_MAX_ARGS_PER_LINE

    for i in xrange(n):
        line = "  "
        for j in xrange(ABAQUS_MAX_ARGS_PER_LINE):
            line += str(dNds.pop()) + ","
        line = line[:-1] + "\n"
        myFile.write(line)
    line = " "
    for i in xrange(r):
        line += str(dNds.pop()) + ","
    line = line[:-1] + "\n\n"
    myFile.write(line)

myFile.write("\n" * 2)

# Frontières
for bId, args1 in boundaries.items():
    for args2 in args1.values():
        myFile.write("*Nset, nset=boundary_" + bId[0] + "_-_" + bId[1] + ", instance=plaque-1\n")

        args2 = args2[0].keys()

        n = len(args2) // ABAQUS_MAX_ARGS_PER_LINE
        r = len(args2) % ABAQUS_MAX_ARGS_PER_LINE

        for i in xrange(n):
            line = "  "
            for j in xrange(ABAQUS_MAX_ARGS_PER_LINE):
                line += str(args2.pop()) + ","
            line = line[:-1] + "\n"
            myFile.write(line)
        line = " "
        for i in xrange(r):
            line += str(args2.pop()) + ","
        line = line[:-1] + "\n\n"
        myFile.write(line)


#for d1, args1 in domainsConnections.items():
    #for d2, args2 in args1.items():
        #myFile.write("*Nset, nset=boundary_" + d1 + "_-_" + d2 + ", instance=plaque_raidie-1\n")

        #args2 = args2.keys()

        #n = len(args2) // ABAQUS_MAX_ARGS_PER_LINE
        #r = len(args2) % ABAQUS_MAX_ARGS_PER_LINE

        #for i in xrange(n):
            #line = "  "
            #for j in xrange(ABAQUS_MAX_ARGS_PER_LINE):
                #line += str(args2.pop()) + ","
            #line = line[:-1] + "\n"
            #myFile.write(line)
        #line = " "
        #for i in xrange(r):
            #line += str(args2.pop()) + ","
        #line = line[:-1] + "\n\n"
        #myFile.write(line)

myFile.write("*End Assembly\n")
myFile.close()






